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Abstract 

Many species show avoidance reactions in response to looming object 
approaches. In locusts, the corresponding escape behavior correlates with 
the activity of the lobula giant movement detector (LGMD) neuron. Dur- 
ing an object approach, its firing rate was reported to gradually increase 
until a peak is reached, and then it declines quickly. The 77-function pre- 
dicts that the LGMD activity is a product between an exponential func- 
tion of angular size exp( — Q) and angular velocity O, and that peak ac- 
tivity is reached before time-to-contact (ttc). The 77-function has become 
the prevailing LGMD model because it reproduces many experimental 
observations, and even experimental evidence for the multiplicative op- 
eration was reported. Several inconsistencies remain unresolved, though. 
Here we address these issues with a new model (-i/j-model) , which explic- 
itly connects and 6) to biophysical quantities. The i/j-model avoids 
biophysical problems associated with implementing exp(-), implements 
the multiplicative operation of 77 via divisive inhibition, and explains why 
activity peaks could occur after ttc. It consistently predicts response fea- 
tures of the LGMD, and provides excellent fits to published experimental 
data, with goodness of fit measures comparable to corresponding fits with 
the 77-function. 



'http://www.ir3c.ub.edu, Research Institute for Brain, Cognition, and Behaviour (IR3C) 
Edifici de Ponent, Campus Mundet, Universitat de Barcelona, Passeig Vail d'Hebron, 171. 
E-08035 Barcelona. Note: A revised version of this paper with the title "Emergence of 
Multiplication in a Biophysical Model of a Wide-Field Visual Neuron for Computing Object 
Approaches: Dynamics, Peaks, & Fits" has been accepted in Advances in Neural Information 
Processing Systems NIPS 2011, Granda, Spain (http://nips.ee). 
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1 Introduction: r and 77 



Collision sensitive neurons were reported in species such different as monkeys 
I], pigeons [361 [34], &ogs [M [2U], and insects |I2 [13 [27J [lOj [38] . This in- 
dicates a high ecological relevance, and raises the question about how neurons 
compute a signal that eventually triggers corresponding movement patterns (e.g. 
escape behavior or interoceptive actions). Here, we will focus on visual stimula- 
tion. Consider, for simplicity, a circular object (diameter 21), which approaches 
the eye at a collision course with constant velocity v. If we do not have any a 
priori knowledge about the object in question (e.g. its typical size or speed), 
then we will be able to access only two information sources. These informa- 
tion sources can be measured at the retina and are called optical variables 
(OVs). The first is the visual angle 0, which can be derived from the number 
of stimulated photoreceptors (spatial contrast). The second is its rate of change 
dQ(t)/dt = ©(£). Angular velocity 8 is related to temporal contrast. 
How should we combine Q and in order to track an imminent collision? The 
perhaps simplest combination is r(t) = 0(i)/0(t) [T3l [18] ■ If the object hit 
us at time t c , then r(t) p» t c — t will give us a running estimation of the time 
that is left until contacfl Moreover, we do not need to know anything about 
the approaching object: The ttc estimation computed by r is independent of 
object size and velocity. Neurons with r-like responses were indeed identified in 
the nucleus retundus of the pigeon brain [34] . In humans, only fast interceptive 
actions seem to rely exclusively on r |371 135] . Accurate ttc estimation, however, 
seems to involve further mechanisms (rate of disparity change |31]). 
Another function of OVs with biological relevance is rj = 0exp(— a©), with a = 
const. 10J . While 77- type neurons were found again in pigeons [34] and bullfrogs 
[20] , most data were gathered from the LGMEo in locusts (e.g. [TOl [01 [71 12"5] ) . 
The r\- function is a phenomenological model for the LGMD, and implies three 
principal hypothesis: (i) An implementation of an exponential function exp(-). 
Exponentation is thought to take place in the LGMD axon, via active membrane 
conductances [8]. Experimental data, though, seem to favor a third-power law 
rather than exp(-). (ii) The LGMD carries out biophysical computations for 
implementing the multiplicative operation. It has been suggested that multi- 
plication is done within the LGMD itself, by subtracting the logarithmically 
encoded variables log© — a© [TpJ[8]. (in) The peak of the 77- function occurs 
before ttc, at visual angle 0(f) = 2arctan(l/a) [9j. It follows ttc for certain 
stimulus configurations (e.g. l/\v\ ;$ 5ms). In principle, t > t c can be accounted 
for by n(t + S) with a fixed delay 5 < (e.g. —27ms). But other researchers 
observed that LGMD activity continuous to rise after ttc even for l/\v\ ^ 5ms 
[2"8] . These discrepancies remain unexplained so far [30J, but stimulation dy- 
namics perhaps plays a role. 

1 This linear approximation gets worse with increasing 0, but turns out to work well until 
short before ttc. 

2 LGMD activity is usually monitored via its postsynaptic neuron, the Descending Con- 
tralateral Movement Detector (DCMD) neuron. This represents no problem as LGMD spikes 
follow DCMD spikes 1:1 under visual stimulation [22] from 300Hz gl] to at least 400Hz [24]. 
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We we will address these three issues by comparing the novel function "ip" with 
the 77-function. 

2 LGMD computations with the ^-function 

A circular object which starts its approach at distance xq and with speed v 
projects a visual angle &(t) = 2 arctan[^/(x — vt)] on the retina [34j [9]. The 
kinematics is hence entirely specified by the half-size-to- velocity ratio l/\v\, and 
x . Furthermore, Q(t) = 2lv/((x - vt) 2 + I 2 ). 

In order to define ip, we consider at first the LGMD neuron as an RC-circuit 
with membrane potential V [T7] 

dV 

C m —=P {V rest - V) + g exc (V exc -V)+ g inh (V mh - V) (1) 
at 

C m = membrane capacitjQ; (3 = 1/R m denotes leakage conductance across the 
cell membrane (R m : membrane resistance); g exc and gi n h are excitatory and 
inhibitory inputs. Each conductance gi (i = rest, exc, ink) can drive the mem- 
brane potential to its associated reversal potential V% (usually Vi n h < V exc ). 
Shunting inhibition means Vi — Vrest- Shunting inhibition lurks "silently" 
because it gets effective only if the neuron is driven away from its resting 
potential. With synaptic input, the neuron decays into its equilibrium state 
Voo = (Vrest/3 + V exc g exc + V lnh g lnh ) / (J3 + g exc + g inh ) according to V(t) = 
Voo(l — exp(—t/T m )). Without external input, V(t 3> 1) — ► V rest - The time 
scale is set by r m . Without synaptic input r m = C m //3. Slowly varying inputs 
9exc,9inh > modify the time scale to approximately r m /(l + (g exc + ginh) IP)- 
For highly dynamic inputs, such as in late phase of the object approach, the 
time scale gets dynamical as well. The i/j-model assigns synaptic inputf^ 

fec(i)=#), Ht) = Ci# - Ai stiro ) + (1 - Ci)e(i) (2a) 

9inh{t) = , = Cotf(f - At stim ) + (1 - Co)O(t) (2b) 

Thus, we say ip(t) = V(t) if and only if g exc and ginh are defined with the last 
equation. The time scale of stimulation is defined by A.t S i m (by default lms). 
The variables and •& are lowpass filtered angular size and rate of expansion, 
respectively. The amount of filtering is defined by memory constants £0 and Ci 
(no filtering for zero). The idea is to continue to generate synaptic input after 
ttc, where <d(t > t c ) — const and thus 0(t > t c ) = 0. Inhibition is first weighted 
by 7, and then potentiated by the exponent e. Hodgkin-Huxley potentiates 
gating variables n,m £ [0, 1] instead (potassium oc n 4 , sodium oc m 3 , [12]) and 
multiplies them with conductances. Gabbiani and co-workers found that the 
function which transforms membrane potential to firing rate is better described 
by a power function with e = 3 than by exp(-) (Figure 4d in [5]). 

3 Here we assume that the membrane potential serves as a predictor for the LGMD's mean 
firing rate. 

4 Set to unity for all simulations 

5 LGMD receives also inhibition from a laterally acting network 1211 . The r\- function con- 
siders only direct feedforward inhibition [22] [B], and so do we. 
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Qe [7.63A°, 160.00A=t temporal resolution At stjn =1.gms l/|v|=20.00ms, 6=1.00,y=7.50, e=3.00,^ =0.90, ^=0.99, n |B|fli =25 




time [ms] time [ms] 



(a) discretized optical variables (b) tp versus 77 

Figure 1: (a) The continuous visual angle of an approaching object is shown 
along with its discretized version. Discretization transforms angular velocity 
from a continuous variable into a series of "spikes" (rescaled) . (b ) The ip function 
with the inputs shown in a, with n re i ax = 25 relaxation time steps. Its peak 
occurs t max = 56ms before ttc (t c = 300ms). An 77 function (a = 3.29) that 
was fitted to ip shows good agreement. For continuous optical variables, the 
peak would occur 4ms earlier, and 77 would have a = AAA with R 2 = 1. For 
n re iax = 10, ip is farther away from its equilibrium at Voo, and its peak moves 
19ms closer to ttc. 



3 Dynamics of the ^-function 

Discretization. In a typical experiment, a monitor is placed a short distance 
away form the insect's eye, and an approaching object is displayed. Computer 
screens have a fixed spatial resolution, and as a consequence size increments of 
the displayed object proceed in discrete jumps. The locust retina is further- 
more composed of a discrete array of ommatidia units. We therefore can expect 
a corresponding step-wise increment of O with time, although optical and neu- 
ronal filtering may smooth O to some extent again, resulting in ■& (figure [1]) . 
Discretization renders discontinuous, what again will be alleviated in For 
simulating the dynamics of ip, we discretized angular size with floor(0), and 
®(t) ps [0(f — Atgtim) — Q(t)]/ At sim . Discretized optical variables (OVs) were 
re-normalized to match the range of original (i.e. continuous) OVs. 
To peak, or not to peak? Rind & Simmons reject the hypothesis that the 
activity peak signals impending collision on grounds of two arguments [28j : (i) 
If 0(t + Atgtim) — 0(t) ~ 3° in consecutively displayed stimulus frames, the 
illusion of an object approach would be lost. Such stimulation would rather be 
perceived as a sequence of rapidly appearing (but static) objects, causing re- 
duced responses, (it) After the last stimulation frame has been displayed (that 
is = const), LGMD responses keep on building up following ttc. This behav- 
ior clearly depends on l/\v\, also according to their own data (e.g. Figure 4 in 
[2"?5]): Response build up beyond ttc is typically observed for sufficiently small 
values of l/\v\. Input into tp in situations where = const and Theta = 0, 
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t c =500ms, dia=1 2.0cm, At 5t|m =1. 00ms, dt=10.00[is, discrete=1 t c -500ms, dia=1 2.0cm, n re|aji =25, dt=10J)0|is, discrete^ 




(a) different n re i ax (b) different At st v 



Figure 2: The figures plot the relative time t max = t c — t of the response peak 
of ip, V(t), as a function of half-size-to-velocity ratio (points). Line fits with 
slope a and intercept S were added (lines). The predicted linear relationship 
in all cases is consistent with experimental evidence ;9:. (a) The stimulus time 
scale is held constant at At s ti m = 1ms, and several LGMD time scales are 
defined by n re i ax . Bigger values of n re i ax move V(t) closer to its equilibrium 
Voc{t), implying higher slopes a in turn, (b) LGMD time scale is fixed at 
n re iax = 25, and At st i m is manipulated. Because of the discretization of optical 
variables (OVs) in our simulation, increasing At st i r n translates to an overall 
smaller number of jumps in OVs, but each with higher amplitude. 



respectively, is accommodated by $ and respectively. 

We simulated (i) by setting At s tim = 5ms, thus producing larger and more 
infrequent jumps in discrete OVs than with At st i m = 1ms (default). As a con- 
sequence, i)(t) grows more slowly (delayed build up of inhibition), and the peak 
occurs later (t rnax = t c — t = 10ms with everything else identical with figure[T]D). 
The peak amplitude V — V(t) decreases nearly sixfold with respect to default. 
Our model thus predicts the reduced responses observed by Rind & Simmons 
[28]. 

Linearity. Time of peak firing rate is linearly related to l/\v\ [TQl [9]. The 
ry-function is consistent with this experimental evidence: t = t c — al/\v\ + 8 
(e.g. a — 4.7, <5 = —27ms). The ^-function reproduces this relationship as well 
(figure [2]), where a depends critically on the time scale of biophysical processes 
in the LGMD. We studied the impact of this time scale by choosing 10/j.s for the 
numerical integration of equation[T] (algorithm: 4th order Runge-Kutta). Apart 
from improving the numerical stability of the integration algorithm, tjj is far 
from its equilibrium Voo(t) in every moment t, given the stimulation time scale 
Atstim — lms[^]. Now, at each value of Q(t) and O(i), respectively, we interca- 
lated n re i ax iterations for integrating ip. Each iteration takes V{t) asymptoti- 
cally closer to Voo(t), and lim nreiax> i V{t) — Voo{t). If the internal processes in 
the LGMD cannot keep up with stimulation (n re i ax =0), we obtain slopes val- 

6 Assuming one At at i m for each integration time step 
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(a) j3 varies 



(b) e varies 



(c) 7 varies 



Figure 3: Each curve shows how the peak ip^ = ipoo(t) depends on the half- 
size-to- velocity ratio. In each display, one parameter of ipoo is varied (legend), 
while the others are held constant (figure title). Line slopes vary according to 
parameter values. Symbol sizes are scaled according to rmse (see also figure 
3]). Rmse was calculated between normalized ipoc(t) & normalized rj(t) (i.e. 
both functions 6 [0, 1] with original minimum and maximum indicated by the 
textbox). To this end, the peak of the ^-function was placed at t c , by choosing, 
at each parameter value, a = \v\ - (t c — (for determining correlation, the mean 
value of a was taken across J/H). 

ues that underestimate experimentally found values (figure [2£i). In contrast, for 
n reiax ~ 25 we get an excellent agreement with the experimentally determined 
a. This means that - under the reported experimental stimulation conditions 
(e.g. [5]) - the LGMD would operate relatively close to its steady statcQ. 
Now we fix n re i ax at 25 and manipulate At st i m instead (figure HJd). The default 
value At st i m = lms corresponds to a — 3.91. Slightly bigger values of At st i m 
(2.5ms and 5ms) underestimate the experimental a. In addition, the line fits 
also return smaller intercept values then. We see t max < up to l/\v\ « 13.5ms 
- LGMD activity peaks after ttc. Or, in other words, LGMD activity continues 
to increase after ttc. In the limit, where stimulus dynamics is extremely fast, 
and LGMD processes are kept far from equilibrium at each instant of the ap- 
proach, a gets very small. As a consequence, t max gets largely independent of 
l/\v\: The activity peak would cling to t max although we varied l/\v\. 

4 Freeze! Experimental data versus steady state 
of "psi" 

In the previous section, experimentally plausible values for a were obtained if ip 
is close to equilibrium at each instant of time during stimulation. In this section 



7 Notice that in this moment we can only make relative statements - we do not have data 
at hand for defining absolute time scales 



(a) j3 varies 



(b) e varies 



(c) 7 varies 



Figure 4: This figure complements figure [3l It visualizes the time averaged 
absolute difference between normalized V'oo(i) & normalized r/(t). For 77, its value 
of a was chosen such that the maxima of both functions coincide. Although not 
being a fit, it gives a rough estimate on how the shape of both curves deviate 
from each other. The maximum possible difference would be one. 




Figure 5: The original data (legend label "HaGaLa95") were resampled from 
ref. [TU] and show DCMD responses to an object approach with = const. 
Thus, increases linearly with time. The ?7-function (fitting function: Ar](t + 5) + 
o) and V'oo (fitting function: Atpooit) + o) were fitted to these data: (a) (Figure 3 
Di in [TU]) Good fits for -0oo are obtained with e — 5 or higher (e = 3 ~> R 2 = 0.35 
and rmse = 0.644; e = 4 ~> R 2 = 0.45 and rmse = 0.592). "Psi" adopts a sigmoid- 
like curve form which (subjectively) appears to fit the original data better than 
rj. (b) (Figure 3 Dii in [10]) "Psi" yields an excellent fit for e = 3. 



we will thus introduce a steady-state version of ip, 

0w + W79(or (3) 
u /3+0(t) + [ 7 0(t)r [> 

(Here we use continuous versions of angular size and rate of expansion). The 
i/'oo-function makes life easier when it comes to fitting experimental data. How- 
ever, it has its limitations, because we brushed the whole dynamic of ip under the 
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Figure 6: (a) DCMD activity in response to a black square {l/\v\ — 30ms, leg- 
end label "e011posl4", ref. [30.) approaching to the eye center of a gregarious 
locust (final visual angle 50°). Data show the first stimulation so habituation is 
minimal. The spike trace (sampled at 10 4 Hz) was full wave rectified, lowpass 
filtered, and sub-sampled to 1ms resolution. Firing rate was estimated with 
Savitzky-Golay filtering ("sgolay"). The fits of the 77-function {Aq{t + 8) + o; 
4 coefficients) and f/'oo-nmction (Aip^t) with fixed e, o, S, Vi n h\ 3 coefficients) pro- 
vide both excellent fits to firing rate, (b) Fitting coefficient a (— > r\- function) 
inversely correlates with f3 (— > ipoo) when fitting firing rates of another 5 trials 
as just described. Similar correlation values would be obtained if e is fixed at 
values e = 2.5,4,5 -~* c — —0.95,-0.96,-0.91. If o was determined by the 
fitting algorithm, then c = —0.70. No clear correlations with a were obtained 
for 7. 



carpet. Figure [3] illustrates how the linear relationship (— "linearity") between 
t max =t c — t and l/\v\ is influenced by changes in parameter values. Changing 
any of the values of e, (3, 7 predominantly causes variation in line slopes. The 
smallest slope changes are obtained by varying Vi n h (data not shown; we checked 
Vinh = 0,-0.001,-0.01,-0.1). For V lnh ^ -0.01, linearity is getting slightly 
compromised, as slope increases with l/\v\ (e.g. V inh — — 1 ~+ a e [4.2,4.7]). 
In order to get a notion about how well the shape of ipoo (t) matches rj(t), we com- 
puted time-averaged difference measures between normalized versions of both 
functions (details: figure [3] & |4| . Bigger values of /3 match 77 better at smaller, 
but worse at bigger values of l/\v\ (figure 0k) . Smaller f3 cause less variation 
across l/\ v \. As to variation of e, overall, curve shapes seem to be best aligned 
with e = 3 to e = 4 (figure 2Jd) . Furthermore, better matches between (t) 
and r]{t) correspond to bigger values of 7 (figure 2h) . And finally, Vinh marches 
again to a different tunc (data not shown). Vinh = — 0.1 leads to the best agree- 
ment (~ 0.04 across l/\v\) of all Vinh-, quite different from the other considered 
values. For the rest, ipoo(t) and r){t) align the same (all have maximum 0.094), 
despite of covering different orders of magnitude with Vi n h = 0, —0.001, —0.01. 
Decelerating approach. Hatsopoulos et al. [10] recorded DCMD activity in 
response to an approaching object which projected image edges on the retina 
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moving at constant velocity: = const, implies Q(t) = &o + Oi. This "linear 
approach" is perceived as if the object is getting increasingly slower. But what 
appears a relatively unnatural movement pattern serves as a test for the func- 
tions 77 & -0OO- Figure [5] illustrates that V'oo passes the test, and consistently 
predicts that activity sharply rises in the initial approach phase, and subse- 
quently declines (77 passed this test already in the year 1995). 
Spike traces. We re-sampled about 30 curves obtained from LGMD recordings 
from a variety of publications, and fitted 77 & ipao -functions. We cannot show 
the results here, but in terms of goodness of fit measures, both functions are 
in the same ballbark. Rather, figure^ shows a representative example [50] . 
When a and (3 are plotted against each other for five trials, we see a strong 
inverse correlation (figure |6j}) . Although five data points are by no means a 
firm statistical sample, the strong correlation could indicate that /3 and a play 
similar roles in both functions. Biophysically, f3 is the leakage conductance, 
which determines the (passive) membrane time constant r m oc 1//3 of the neu- 
ron. Voltage drops within r m to its 1/e part. Bigger values of /3 mean shorter 
r m (i.e., "faster neurons"). Getting back to 77, this would suggest a oc r m , such 
that higher (absolute) values for a would possibly indicate a slower dynamic of 
the underlying processes. 

5 Discussion ("The Good, the Bad, and the Ugly") 

Up to now, mainly two classes of LGMD models existed: The phenomenological 
77-function on the one hand, and computational models with neuronal layers 
presynaptic to the LGMD on the other (e.g. [351 [TS]; real- world video sequences 

6 robotics: e.g. [3J Q31 [351 [2]). Computational models predict that LGMD 
response features originate from excitatory and inhibitory interactions in - and 
between - presynaptic neuronal layers. Put differently, non-linear operations are 
generated in the presynaptic network, and can be a function of many (model) 
parameters (e.g. synaptic weights, time constants, etc.). In contrast, the 77- 
function assigns concrete nonlinear operations to the LGMD [7] . The 77-function 
is accessible to mathematical analysis, whereas computational models have to 
be probed with videos or artificial stimulus sequences. The 77-function is vague 
about biophysical parameters, whereas (good) computational models need to 
be precise at each (model) parameter value. The 77-function establishes a clear 
link between physical stimulus attributes and LGMD activity: It postulates 
what is to be computed from the optical variables (OVs). But in computational 
models, such a clear understanding of LGMD inputs cannot always be expected: 
Presynaptic processing may strongly transform OVs. 

The ip function thus represents an intermediate model class: It takes OVs as 
input, and connects them with biophysical parameters of the LGMD. For the 
neurophysiologist, the situation could hardly be any better. Psi implements the 
multiplicative operation of the 77-function by shunting inhibition (equation [1} 
Vexc ~ V res t and Vi n h ~ V res t). The 77-function fits ip very well according to our 
dynamical simulations (figure [lj, and satisfactory by the approximate criterion 
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of figure |H 

We can conclude that tp implements the 77-function in biophysically plausible 
way. However, ip does neither explicitly specify rj's multiplicative operation, 
nor its exponential function exp(-). Instead we have an interaction between 
shunting inhibition and a power law (-) e , with e 3. So what about power laws 
in neurons? 

Because of e > 1, we have an expansive nonlinearity. Expansive power- law 
nonlinearities are well established in phenomenological models of simple cells 
of the primate visual cortex [TJ [TT] . Such models approximate a simple cell's 
instantaneous firing rate r from linear filtering of a stimulus (say V) by r oc 
([Y] + ) e , where [•]+ sets all negative values to zero and lets all positive pass. 
Although experimental evidence favors linear thresholding operations like r cx 
[Y — Y t hres} + , neuronal responses can behave according to power law functions 
if Y includes stimulus- independent noise [15]. Given this evidence, the power- 
law function of the inhibitory input into ip could possibly be interpreted as a 
phenomenological description of presynaptic processes [2] . 
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